###Figure of chlorophyll absorption

setwd('/media/Iomega_HDD/boulot/Article/kd/last_23_01_2012/data/')
wave<-c(320,340,380,412,443,490,510,555)
nwave<-length(wave)
kw_MM3<-read.table('asw_bbsw_300_1100nm.txt',header=T,sep="\t")

kw <-kw_MM3[which(kw_MM3[,1]%in%wave),2]+kw_MM3[which(kw_MM3[,1]%in%wave),3]

load("dataset.Rdata")
table <- cbind(signif(t5$value,7), signif(t5$gsm_chl,7), signif(t5$gsm_cdm,7),
                   signif(t5$gsm_bbp,7),t5$lambda, t5$convert_lambda)

para<-read.table('../script_admb/model.std',header=T, fill=T)
para <- para[,1:4]

s0<-para[9,3]
aphstar<-para[27:34,3]
a2<-para[18:25,3]
eta<-para[26,3]


ch25<-quantile(table[which(table[,6]==6),2],0.25)
ch50<-quantile(table[which(table[,6]==6),2],0.5)
ch75<-quantile(table[which(table[,6]==6),2],0.75)
meanch<-mean(table[which(table[,6]==6),2])

acdm25<-quantile(table[which(table[,6]==6),3],0.25)
acdm50<-quantile(table[which(table[,6]==6),3],0.5)
acdm75<-quantile(table[which(table[,6]==6),3],0.75)
meanacdm<-mean(table[which(table[,6]==6),3])
acd<-c(round(acdm25,digits=3),round(acdm50,digits=2),round(acdm75,digits=2))#,round(meanacdm,digits=2))

bb25<-quantile(table[which(table[,6]==6),4],0.25)
bb50<-quantile(table[which(table[,6]==6),4],0.5)
bb75<-quantile(table[which(table[,6]==6),4],0.75)
meanbb<-mean(table[which(table[,6]==4),3])
bbb<-c(round(bb25,digits=4),round(bb50,digits=4),round(bb75,digits=4))#,round(meanbb,digits=4))

################
abric<-exp(c(NA,NA,NA,log(0.0323),log(0.039),log(0.0274),log(0.018),log(0.0068)))
a2bric<-c(NA,NA,NA,0.714,0.652,0.632,0.74, 0.9732)


a<-c(round(ch25,digits=2),round(ch50,digits=1),round(ch75,digits=1),round(meanch,digits=1),3,10)
a3<-matrix(0,nrow=length(a),ncol=nwave)
a4<-matrix(0,nrow=length(a),ncol=nwave)
for (j in 1:length(a)){
for (i in 1:nwave){
a3[j,i]<-aphstar[i]*a[j]^(a2[i])
a4[j,i]<-abric[i]*a[j]^a2bric[i]
}}

opar<-par()
postscript('../Figures/Figure4.ps', height=10, width=14, paper="special", horizontal=F)
color<- terrain.colors(length(a)+6)
color2<- heat.colors(length(a)+6)
par(mar=c(5,6,1,1))
par(mfrow=c(2,2))
plot(a3[1,]~wave, type="o", ylim=c(0,max(a3)), col=color[1], lwd=3, pch=16, lty=1, cex.axis=2, cex.lab=2, ylab=expression(paste("Chlorophyll absorption", "  (",m^-1,")")), xlab=expression(paste(lambda, "(nm)")), cex=1.2)
text(550, 0.35, "(a)", cex=2)
lines(a4[1,]~wave, type="o", ylim=c(0,max(a3)), col=color2[1], lwd=3, pch=16, lty=2, cex=1.2)
#plot(a4[1,]~wave, type="o", ylim=c(0,max(a3)), col=color[1], lwd=2, pch=16, lty=2, cex=1.2)
for (i in 2:length(a)){
lines(a3[i,]~wave, type="o", col=color[i], lwd=3, pch=(15+i), lty=1, cex=1.2)
lines(a4[i,]~wave, type="o", col=color2[i], lwd=3, pch=(15+i), lty=2, cex=1.2)
#lines(a4[i,]~wave, type="o", col=color[i], lwd=2, pch=(15+i), lty=2, cex=1.2)
}
legend(310.5,0.36, col=color2[1:length(a)], lwd=rep(2,length(a)), legend=c("","","","","",""), bty="n", cex=2, pch=seq(16,(15+length(a)),by=1), lty=seq(1,length(a)))
legend("topleft", col=color[1:length(a)], lwd=rep(2,length(a)), legend=c(paste("Chl= ",round(ch25,digits=2)),paste("Chl= ",round(ch50,digits=1)),paste("Chl= ",round(ch75,digits=1)),paste("Chl= ",round(meanch,digits=1)),"Chl= 3","Chl= 10"), bty="n", cex=2, pch=seq(16,(15+length(a)),by=1), lty=seq(1,length(a)))


plot((acd[length(acd)]*exp(-s0*(wave-443)))~ wave, col="orange2", lwd=3, pch=16, lty=1, cex.axis=2, cex.lab=2, ylab=expression(paste(a[CDM], "  (",m^-1,")")), xlab=expression(paste(lambda, "(nm)")), cex=2)
text(350, 0.27, "(b)", cex=2)
curve(acd[length(acd)]*exp(-s0*(x-443)),310,560,col="orange2", lwd=3,lty=1,ylab="", xlab="", add=T)

points((acd[length(acd)-1]*exp(-s0*(wave-443)))~ wave, col="orange", lwd=3, pch=17, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
curve(acd[length(acd)-1]*exp(-s0*(x-443)),310,560,col="orange", lwd=3,lty=2,ylab="", xlab="", add=T)

points((acd[1]*exp(-s0*(wave-443)))~ wave, col="gold", lwd=3, pch=18, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
curve(acd[1]*exp(-s0*(x-443)),310,560,col="gold", lwd=3,lty=3,ylab="", xlab="", add=T)

legend("topright", col=c("orange2", "orange", "gold"), lwd=rep(3,length(acd)), legend=c(expression(paste(a[CDM](443),"= 0.009")), expression(paste(a[CDM](443),"= 0.02")),expression(paste(a[CDM](443),"= 0.05"))), bty="n", cex=2, pch=c(16,17,18), lty=seq(1,3))

########For bbp
plot((bbb[3]*(wave/443)^eta)~ wave, col="grey68", lwd=3, pch=16, lty=1, cex.axis=2, cex.lab=2, ylab=expression(paste(b[bp], "  (",m^-1,")")), xlab=expression(paste(lambda, "(nm)")), ylim=c(0,0.005), cex=2)
text(350, 0.00485, "(c)", cex=2)
curve(bbb[3]*(x/443)^eta,310,560,col="grey68", lwd=3,lty=1,ylab="", xlab="", add=T)

points((bbb[2]*(wave/443)^eta)~ wave, col="grey45", lwd=3, pch=17, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
curve(bbb[2]*(x/443)^eta,310,560,col="grey45", lwd=3,lty=2,ylab="", xlab="", add=T)

points((bbb[1]*(wave/443)^eta)~ wave, col="grey22", lwd=3, pch=18, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
curve(bbb[1]*(x/443)^eta,310,560,col="grey22", lwd=3,lty=3,ylab="", xlab="", add=T)

legend("topright", col=c("grey22", "grey45", "grey68"), lwd=rep(3,length(acd)), legend=c(expression(paste(b[bp](443),"= 0.0011")), expression(paste(b[bp](443),"= 0.0015")),expression(paste(b[bp](443),"= 0.003"))), bty="n", cex=2, pch=c(16,17,18), lty=seq(1,3))

########For Kd
kd1<-(bbb[1]*(wave/443)^eta+kw+acd[1]*exp(-s0*(wave-443))+a3[1,])
kd2<-(bbb[2]*(wave/443)^eta+kw+acd[2]*exp(-s0*(wave-443))+a3[2,])
kd3<-(bbb[3]*(wave/443)^eta+kw+acd[3]*exp(-s0*(wave-443))+a3[3,])


plot(kd1~ wave, col="grey22", lwd=3, pch=16, lty=1, cex.axis=2, cex.lab=2, ylab=expression(paste(K[d](lambda), " (",m^-1,")")), xlab=expression(paste(lambda, "(nm)")), ylim=c(0,0.4), cex=2)
text(350, 0.385, "(d)", cex=2)
tutu<-spline(wave,kd1, n=100)
lines(tutu$y~tutu$x, lwd=3, lty=1, col="grey22")

points(kd2~ wave, col="grey45", lwd=3, pch=17, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
tutu<-spline(wave,kd2, n=100)
lines(tutu$y~tutu$x, lwd=3, lty=1, col="grey45")

points(kd3~ wave, col="grey68", lwd=3, pch=18, cex.axis=2, cex.lab=2, ylab="", xlab="", cex=2)
tutu<-spline(wave,kd3, n=100)
lines(tutu$y~tutu$x, lwd=3, lty=1, col="grey68")

legend("topright", col=c("grey22", "grey45", "grey68"), lwd=rep(3,length(acd)), legend=c(expression(paste(1^st," Quartile")), "Median",expression(paste(3^rd," Quartile"))), bty="n", cex=2, pch=c(16,17,18), lty=seq(1,3))



par(opar)
dev.off()
